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Abstract 

Background: The swimming crab, Portunus trituberculatus, which is naturally distributed in the coastal waters of 
Asia-Pacific countries, is an important farmed species in China. Salinity is one of the most important abiotic factors 
that influence not only the distribution and abundance of crustaceans, it is also an important factor for artificial 
propagation of the crab. To better understand the interaction between salinity stress and osmoregulation, we 
performed a transcriptome analysis in the gills of Portunus trituberculatus challenged with salinity stress, using the 
lllumina Deep Sequencing technology. 

Results: We obtained 27,696,835, 28,268,353 and 33,901,271 qualified lllumina read pairs from low salinity 
challenged (LC), non-challenged (NC), and high salinity challenged (HC) Portunus trituberculatus cDNA libraries, 
respectively. The overall de novo assembly of cDNA sequence data generated 94,511 unigenes, with an average 
length of 644 bp. Comparative genomic analysis revealed that 1 ,705 genes differentially expressed in salinity stress 
compared to the controls, including 615 and 1,516 unigenes in NC vs LC and NC vs HC respectively. GO functional 
enrichment analysis results showed some differentially expressed genes were involved in crucial processes related 
to osmoregulation, such as ion transport processes, amino acid metabolism and synthesis processes, proteolysis 
process and chitin metabolic process. 

Conclusion: This work represents the first report of the utilization of the next generation sequencing techniques for 
transcriptome analysis in Portunus trituberculatus and provides valuable information on salinity adaptation 
mechanism. Results reveal a substantial number of genes modified by salinity stress and a few important salinity 
acclimation pathways, which will serve as an invaluable resource for revealing the molecular basis of osmoregulation 
in Portunus trituberculatus. In addition, the most comprehensive sequences of transcripts reported in this study 
provide a rich source for identification of novel genes in the crab. 
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Introduction 

Portunus trituberculatus (Crustacea: Decapoda: Brachyura), 
commonly known as the swimming crab, is widely distributed in 
the coastal waters of Korea, Japan, China, and southeast Asia 
[1]. This species inhabits estuaries and coastal waters, which 
belong to typical euryhaline crab species. In China, it is a major 
edible crab species and one of the most important fishery 
resources [2] and the production has now reached 90,000 tons, 
valued at AUS$2.5 billion in 2009 [3]. 



salinity is one of the most important abiotic factors that 
influence not only the distribution and abundance of 
crustaceans, but also their general physiology and well being 
[4]. The water salinity condition is also an important factor for 
artificial propagation of the swimming crab [5]. Throughout their 
prolonged culture period, Portunus trituberculatus often 
experience substantial salinity fluctuations either due to heavy 
rainfalls or droughts, which could have significant impacts to 
farm productivity and in severe situations, heavy mortality. This 
often requires the ability to regulate hemolymph osmolytes in 
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relation to the environment they inhabit via osmoregulation to 
control their hemolymph osmotic pressure [6]. 

Due to the implications and critical importance of 
osmoregulation to the crab artificial propagation, a number of 
researchers have been devoted to this topic. An extensive 
literature that describes the growth, development, physiology, 
behavior, and propagation techniques of Portunus 
trituberculatus exposed to salinity stress have revealed the 
crab grows in optimal salinity ranged from 20-35ppt, whereas 
they can occur at salinities below 6 ppt and will survive 
salinities in excess of 48ppt [7-12]. In order to study the 
mechanism of osmoregulation, Xu et.al. investigated gene 
expression in the Portunus trituberculatus exposed to different 
salinity stresses via cDNA microarray chip, and 417 
differentially expressed genes were identified [5]. Their study 
revealed a few important salinity acclimation pathways, which 
may be helpful in understanding the molecular basis of 
osmoregulation and salinity adaptation in the crab. 

Even though cDNA microarray technology is a powerful tool 
for studying genome-wide gene expression, this technology 
fails to detect sequence variation and to recognize new genes 
or transcripts and can only be designed from limited expressed 
sequence tag data as the genome of Portunus trituberculatus 
has not yet been determined. To date, there are 13,985 ESTs 
available for the crab in the Genebank, however, it remains 
insufficient for the comprehensive understanding of Portunus 
trituberculatus transcriptome. Many low expression transcripts 
would be missed from current EST data, which makes it difficult 
for further analysis on transcriptome. 

Next-generation high-throughput RNA sequencing 
technology (RNA-seq) is a recently-developed method for 
discovering, profiling, and quantifying RNA transcripts with 
several advantages over other expression profiling 
technologies including higher sensitivity and the ability to detect 
splicing isoforms and somatic mutations [13]. Because it is not 
restricted by the unavailability of a genome reference 
sequence, this approach has been applied in decoding the 
genomes of several non-model organisms, providing valuable 
information in the understanding of gene function, cell 
responses and evolution [14-16]. Significant progress has also 
been made in understanding the transcript expression of 
various marine Crustacea by RNA-seq over the last two years, 
such as Litopenaeus vannamei, Fenneropenaeus chinensis, 
Eriocheir sinensis and Macrobrachium nipponense [17-21]. The 
countable, almost digital, nature of RNA-seq data makes them 
particularly attractive for the quantitative analysis of transcript 
expression levels, which can give reliable measurements of 
transcript levels in one or more conditions [22]. However, such 
investigations in Portunus trituberculatus have not been 
reported. 

In the present study, we examined the whole transcriptome 
responses to salinity stress of the Portunus trituberculatus for 
the first time using the lllumina's sequencing technology. 
Considering individual monitoring of the Portunus 
trituberculatus responses to salinity stress, nine libraries (three 
technical replicates of each condition) were established from 
the gills of Portunus trituberculatus that exposed to optimal, low 
and high salinity seawater, respectively. The study aimed to 



compare the expression patterns of the three conditions to 
better understand the transcriptomic regulation in Portunus 
trituberculatus to salinity stress and identify genes involved in 
osmoregulation of the crab. The results of this study are an 
important resource for future researches on mechanism of 
osmoregulation for marine invertebrates. 

Materials and Methods 

Salinity Challenge Experiment and Sample Preparation 

The swimming crabs, Portunus trituberculatus at 80 days 
age (5.62-1 1 .66g in body weight), were obtained from a local 
farm in Qingdao, China. All the samples were acclimated in the 
laboratory (33 ppt, 18°C) for one week before the experiment 
treatment. The crabs were divided into 3 groups (90 crabs for 
one group) and acclimated to low salinity challenged (LC, 5 
ppt), non-challenged (NC, 33 ppt), high salinity challenged (HC, 
50 ppt) at 18°C. For each group, the sixth pair of gills from 9 
crabs were collected after ten days and samples were stored in 
RNAIater (Ambion) at 4 °C over night and then at -20 °C until 
RNA extraction within 2 weeks. 

RNA Isolation, cDNA Library Construction and lllumina 
Deep Sequencing 

Total RNA was isolated from each sample by trizol 
(Invitrogen, CA,USA). RNA degradation and contamination was 
monitored on 1% agarose gels. RNA purity was checked using 
the Nano Photometer® spectrophotometer (IMPLEN, CA, 
USA) . RNA concentration was measured using Qubit® RNA 
Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, 
USA). RNA integrity was assessed using the RNA Nano 6000 
Assay Kit of the Bioanalyzer 2100 system (Agilent 
Technologies, CA, USA). A total amount of 3g RNA per sample 
was used as input material for the RNA sample preparations 
and all samples had RIN values above 8. Then, samples of 
three individuals were pooled within each group in equal 
amounts to generate three mixed sample. The pooling samples 
were then used to prepare 9 separate lllumina sequencing 
libraries (three technical replicates of each condition). 

cDNA libraries were generated using lllumina TruSeq™ RNA 
Sample Preparation Kit (lllumia, San Diego, USA) following 
manufacturer's recommendations and index codes were added 
to attribute sequences to each sample. Briefly, mRNA was 
purified from total RNA using poly-T oligo-attached magnetic 
beads. Fragmentation was carried out using divalent cations 
under elevated temperature in lllumina proprietary 
fragmentation buffer. First strand cDNA was synthesized using 
random oligonucleotides and Superscript II. Second strand 
cDNA synthesis was subsequently performed using DNA 
Polymerase I and RNase H. Remaining overhangs were 
converted into blunt ends via exonuclease/polymerase 
activities and enzymes were removed. After adenylation of 3' 
ends of DNA fragments, lllumina PE adapter oligonucleotides 
were ligated to prepare for hybridization. In order to select 
cDNA fragments of preferentially 200 bp in length the library 
fragments were purified with AMPure XP system (Beckman 
Coulter, Beverly, USA). DNA fragments with ligated adaptor 
molecules on both ends were selectively enriched using 
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lllumina PCR Primer Cocktail in a 10 cycle PCR reaction. 
Products were purified (AMPure XP system) and quantified 
using the Agilent high sensitivity DNA assay on the Agilent 
Bioanalyzer 2100 system. 

In the final step before sequencing, all 9 individual libraries 
were normalised and pooled together in a single lane on an 
lllumina HiSeq2000 platform and 90 bp paired-end reads were 
generated. 

Bioinformatic Analysis 

Quality control. Raw data (raw reads) of fastq format were 
firstly processed through our self-written perl scripts. In this 
step, clean data (clean reads) were obtained by removing 
reads containing adapter, reads containing ploy-N and low 
quality reads from raw data. At the same time, Q20, Q30, GC- 
content and sequence duplication level of the clean data were 
calculated. All the downstream analyses were based on clean 
data with high quality. 

Transcriptome assembly. The left files (readl files) from 
all libraries/samples were pooled into one big left.fq file, and 
right files (read2 files) into one big right.fq file. Transcriptome 
assembly was accomplished based on the left.fq and right.fq 
using Trinity [23] with min_kmer_cov set to 2 and all other 
parameters set default. 

Transcriptome Annotation. The unigenes were compared 
with the protein nonredundant database using BlastX with E 
values less than 1.0x10 5 (E values less than 1.0x10 5 were 
considered as significant) [24]. All annotated unigenes were 
used to determine the Clusters of Orthologous Groups of 
proteins (COG) term, Gene Ontology (GO) term and the Kyoto 
Encyclopedia of Genes and Genomes (KEGG) pathway with a 
cut-off E-value of 1.0x1 0* 5 using BlastX. 

Differential Expression, Cluster analysis and GO 
enrichment Analysis. Differential expression analysis was 
performed using the DEGseq (2010) R package. P value was 
adjusted using q value [25]. qvalue<0.005&|log2 (foldchange)| 
>1was set as the threshold for significantly differential 
expression. The three samples from each treatment were used 
to generate mean expression levels. Hierarchical cluster 
analysis is used to identify differentially expressed genes with 
certain patterns of expression under two different salinity 
challenge using Cluster 3.0 [5,26]. The settings for the 
calculations were as follows [5]: similarity is measured by 
standard correlation, and clustering method is average linkage. 
The genes are clustered according to similarities in expression 
profiles across three salinity conditions. GO enrichment 
analysis of the differentially expressed genes (DEGs) was 
performed using GOseq based Wallenius non-central hyper- 
geometric distribution[27], which can adjust for gene length 
bias in DEGs. 

Real-time RT-PCR confirmation of lllumina sequencing 
data 

To validate our lllumina sequencing data, twelve differentially 
expressed genes were selected for quantitative RT-PCR 
analysis, using the same RNA samples as for the 
transcriptome profiling. Primers were designed using the 
Primer5 software (Premier Biosoft International) (Table S1). 



RpL8 gene was selected as an internal control for qPCR 
analysis and the primers reference Xu's literature [5]. First 
strand cDNA was synthesized froml mg of RNA using M-MuLV 
reverse transcriptase (Qiagen). The qPCR reaction mixture (20 
mL) consisted of 26 Power SYBR Green PCR Master mix, 
0.9M each of the forward and reverse primers, and 1 mL of 
template cDNA. PCR amplification was performed under the 
following conditions: 50°C for 2 min and 95°C for 30 s, followed 
by 40 cycles of 95°C for 15 s and 62°C for 1 min, and a final 
extension at 72°C for 5 min. 

Result 

lllumina Draft Reads and Sequence Assembly 

In this study, gills of three crabs were used to prepare one 
pooled RNA sample for each group of LC, NC and HC (three 
technical replicates of each group). Nine cDNA libraries were 
then constructed to perform lllumina sequencing. The 
schematic of lllumina deep sequencing and analysis are shown 
in Figure 1. We obtained 27,696,835, 28,268,353, and 
33,901,271 qualified lllumina read pairs for LC, NC and HC, 
giving rise to total clean bases of 5.60, 5.72 and 6.80G, 
respectively. The overall lllumina read pairs and clean bases 
for all samples are 89,866,459 and 18.16G, respectively (Table 
1). Files containing these data were deposited in the Short 
Read Archive of the National Center for Biotechnology 
Information (NCBI) with accession numbers of SRR1013694 
(LC), SRR1013695 (NC) and SRR1013696 (HC). 

After assembly analysis based on all lllumina reads, we 
identified 141,339 contigs. The average length of all contigs 
was 1,157 bp, with the smallest sequence of 201 bp and the 
largest one of 36,385 bp. The sequence length distribution of 
contigs is indicated in Figure 2 and Table 2. Further assembly 
analysis showed that all contigs contributed to 94,511 
unigenes, with an average length of 644 bp. 

To assess the abundance and coverage of the transcriptome 
data, we matched the assembled unigenes against the known 
EST library on Genbank. The 13,985 ESTs downloaded from 
NCBI were clustered and assembled, and 2,612 assembled 
EST-unigenes with mean length of 783 bp were generated. 
Comparisons between transcriptome unigenes and EST- 
unigenes were performed using BLASTn algorithm. Results 
were shown in Figure 3 as a Venn chart. 74.3% (1,942 out of 
2,612) of the EST-unigenes can be matched in the 
transcriptome unigenes library, whereas only 2.05% (1 ,942 out 
of 94,511) of the transcriptome unigene sequences can be 
found in the ESTs library. 

Annotation of Unigenes 

After ruling out short-length and low-quality sequences, 
94,511 unigenes were selected and subjected to annotation 
analysis by matching sequences against Nr and Swiss-prot 
databases using BLASTx searching with an E value 1.0x10 5 
and all annotated unigenes were used to determine the COG 
term, GO term and KEGG pathway (Table 3 and Figure S1, S2, 
S3). 13,212 unigenes (13.98% of the total) can be matched in 
Nr database, and 11,711 (12.39% of the total) matched in 
Swissprot. For main species distribution matched against Nr 
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Figure 1. Schematic of Illumina deep sequencing and analysis. It includes sample preparation, cDNA library construction and 
Illumina sequencing, data analysis including assemble, blast, GO annotation, gene expression analysis, etc. 

doi: 10.1371/joumal.pone.0082155.g001 
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Figure 3. Venn chart for comparisons between assembled transcriptome unigenes and assembled EST-unigenes. 

doi: 10.1371/joumal.pone.0082155.g003 



Table 1. Summary of draft reads of three libraries by 
lllumina deep sequencing. 



Group a 


PE library 


size (bp) Read pairs 


Read length (bp) 


Cleanbases (G) 


LC 


200 


27,696,835 


90 


5.60 


NC 


200 


28,268,353 


90 


5.72 


HC 


200 


33,901,271 


90 


6.84 


Total 


200 


89,866,459 


90 


18.16 



a, the data of every group comes from three technical replicates 
doi: 10.1371/joumal.pone.0082155.t001 



Table 2. Assembly Results Statistics. 





Transcript length interval 


200-500bp 


500-1 kbp 1k-2kbp >2kbp 


total 


Number of transcripts 


73,467 


24,848 19,568 23,456 


141,339 



doi: 10.1371/journal.pone.0082155.t002 

database, 1 1 .25% of the matched unigenes showed similarities 
with Daphnia pulex, followed by Tribolium castaneum (6.90%), 
Pediculus humanuscorporis (5.41%), Branchiostoma floridae 
(3.98%), Paramecium tetraurelia (3.58%), Strongylocentrotus 



purpuratus (3.42%), Nasonia vitripennis (3.17%), Ixodes 
scapularis (2.90%), Megachile rotundata (2.84%), 
Acyrthosiphon pisum (2.66%), Camponotus floridanus (2.1 1 %), 
Harpegnathos saltator (2.07%), Danaus plexippus (1.98%), 
Pediculus humanus (1.57%), Tetrahymena thermophila 
(1.36%), Aedes aegypti (1.36%), and others (43.45%). 

Differentially Expressed Genes 

Comparison of gene expression showed that a total of 1 ,705 
unigenes were differentially expressed after two different salt 
challenges (qvalue<0.005 & |log2 (foldchange)|>1), including 
615 (158 up and 457 down) differentially expressed unigenes 
between NC and LC (Table S2), and 1516 (895 up and 621 
down) between NC and HC (Table S3). Moreover, 426 
unigenes were significantly differentially expressed in both NC 
Vs. LC and NC Vs. HC (Figure 4). To validate our lllumina 
sequencing results, tweleve differentially regulated genes were 
selected for quantitative real time-PCR (qRT-PCR) analysis, of 
which ten genes agrees well with the illumina sequencing 
analysis. The other two, including sodium-potassium-chloride 
cotransporter and V-type proton ATPase subunit F, do not 
match the lllumian sequencing analysis perfectly, however, it 
shows similar trends in up- or down-regulation of genes 
analysised by lllumina sequencing (table 4). The results 
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Table 3. Summary statistics of Portunus trituberculatus transcriptome blast assignment. 





Unigenes 


Nrannotations 


Swiss-protannotations 


ESTscan prediction 


COGhits 


GOmapped 


KEGG hits 


94,511 


13,211 


11,711 


39,625 


1 1 ,528 


23,914 


5,419 



doi: 10.1371/journal.pone.0082155.t003 



LC vs NC 



HC vs NC 




Figure 4. Differentially expressed genes that are unique or shared among LC vs NC and HC vs NC. 

doi: 10.1371/journal.pone.0082155.g004 



suggesting that the majority of the salinity responsive genes 
found from the lllumina sequencing data analysis are authentic. 

The differentially expressed genes were further categorized 
into eight patterns based on expression profiles via hierarchical 
cluster analysis (Figure 5): genes highly up-regulated in low 
and high salinity conditions (Cluster I, 39 genes), only in low 
salinity challenge (Cluster V, 104 genes) and only in high 
salinity conditions (Cluster VII, 840 genes). These categories 
are in contrast with three clusters that showed different kinetics 
of down-regulation in the two salinity conditions (Cluster II, 356 
genes), only in low salinity challenge (Cluster VI, 85 genes), 
and only in high-salinity environment (Cluster VIII, 250 genes). 
Moreover, Cluster III (15 genes) showed the genes highly up- 
regulated in low salinity and down-regulated in high salinity 
challenge, while Cluster IV (16 genes) exhibited the genes 



highly up-regulated in high salinity and down-regulated in low 
salinity challenge. 

GO enrichment analysis 

Biological process of the differentially expressed genes of 
the eight patterns were further analyzed according to GO 
functional enrichment analysis (Figure 6, Table 5 and Table 
S4). GO functional enrichment analysis results showed cluster I 
were involved in chitin metabolic process (GO:0006030, 3 
unigenes). Cluster II were involved in ten biological processes 
at a significant level. The mainly processes and involved 
unigene numbers are chitin metabolic process (GO:0006030, 
13 unigenes), protein polymerization (GO:0051258, 6 
unigenes), response to oxidative stress (GO:0006979, 4 
unigenes), potassium ion transport (GO:0006813, 4 unigenes). 
Cluster III were involved in four biological processes which 
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Table 4. Gene expression levels of 12 selected genes for 
the comparison of controland salinity-challenged groups 
(5% or 45%)according to lllumina sequence and qPCR 
analysis. 













lllumina 




Gene 


ID 


qPCR 




sequence 






LC vs. 


HC vs. 


LC vs. 


HC vs. 






NC 


NC 


NC 


NC 


glutamine 
synthetase 


comp33831 cO 


0.81 


0.20* 


0 73 


0 16* 


glutamate 
dehydrogenase 


comp40590 cO 


1.10 


0.31* 


1 20 


0 36* 


Sodium / potassium 
ATPase beta chain 


comp53011 cO 


0.1 1* 


0.27* 


0 13* 


0 25* 


aquaporin 


comp54248_c0 


0.32* 


0.10* 


0.23* 


0.01* 


heat shock protein 
70 


comp54992_c0 


0.53 


4.26* 


0.64 


2.35* 


carbonic anhydrase 


comp55558_c0 


7.36* 


0.54* 


4.73* 


0.44* 


V-type proton 
ATPase subunit B 


comp62089_c0 


1.11 


0.39* 


0.59 


0.46* 


sodium-potassium- 












chloride 


comp70314_c0 


4.13* 


0.31* 


1.03 


0.26* 


cotrans porter 












chitinase 


comp71131_c0 


0.34* 


0.16* 


0.23* 


0.04* 


V-type proton 
ATPase subunit F 


comp74516_c3 


3.46* 


5.39* 


1.45 


2.13* 


Chloride channel 
protein 


comp76355_c5 


6.71* 


1.63 


2.80* 


0.85 


sodium/hydrogen 
exchanger 


comp76376_c2 


9.10* 


1.35* 


2.56* 


1.81 



Asterisk (*) marks the significant difference betweensalinity-challenged group and 

control group (P < 0.05). 

doi: 10.1371/journal.pone.0082155.t004 



mainly include oxidation-reduction process (GO:0055114, 4 
unigenes), lipid metabolic process (GO:0006629, 3 unigenes). 
Cluster IV were involved in five biological processes which 
mainly include oxidation-reduction process (GO:0055114, 5 
unigenes), proteolysis (GO:0006508, 3 unigenes). Cluster V 
were involved in three biological processes which mainly 
include defense response (GO:0006952, 5 unigenes), anion 
transport (GO:0006820, 5 unigenes). Cluster VI were involved 
in five biological processes which mainly include proteolysis 
(GO:0006508, 8 unigenes), chitin metabolic process (GO: 
0006030, 3 unigenes). Additionally, cluster VII and cluster VIII 
were involved in twenty-nine and sixteen biological processes 
at a significant level respectively. The mainly processes, 
involved unigene numbers and cluster numbers are oxidation- 
reduction process (GO:0055114, 72 unigenes, cluster VII), 
cellular protein modification process (GO:0006464, 50 
unigenes, cluster VII), proteolysis (GO:0006508, 40 unigenes, 
cluster VII), phosphorylation (GO:0016310, 31 unigenes, 
cluster VII), oxidation-reduction process (GO:0055114, 27 
unigenes, cluster VIII), proteolysis (GO:0006508, 20 unigenes, 
cluster VIII), carbohydrate metabolic process (GO:0005975, 17 



unigenes, cluster VIII), chitin metabolic process (GO:0006030, 
10 unigenes, cluster VIII). 

Discussion 

In recent years, the next generation sequencing methods 
have also been applied to analyze transcriptomes of 
crustaceans using lllumina sequencing. Comparing with the 
traditional methods, the next-generation high-throughput DNA 
sequencing techniques provide more ideal methods for 
transcriptome analyses with high efficiency, low cost and high 
data output. The former studies on Portunus trituberculatus 
transcriptome were performed using traditional cDNA library 
and Sanger sequencing methods with RNA from many organs 
such as gill, eyestalk, blood and hepatopancreas [28-31], 
however, it remains insufficient for the comprehensive 
understanding of Portunus trituberculatus transcriptome. In our 
study, we generated 141,339 contigs of Portunus 
trituberculatus transcriptome based on the next generation 
sequencing techniques. Further assembly analysis showed that 
all contigs contributed to 94,511 unigenes, largely enriching the 
transcriptome data of Portunus trituberculatus and prompting 
the genome studies of crustaceans. 

We compared our transcriptome data with Portunus 
trituberculatus EST sequences obtained from NCBI and 
showed that more than half of the EST sequences (74.3%) can 
be matched in the transcriptome data, whereas only 2.05% of 
the transcriptome unigene sequences can be found in the 
ESTs library. It suggests the transcriptome data provide 
abundant information besides the now available ESTs 
sequences. After a homology search in the non-redundant 
protein database at NCBI, a total of 13,212 unigenes, which 
took up a proportion of 13.98% in all the unigenes, showed 
significant BlastX hits of known protein sequences. The 
distribution of significant BlastX hits over different organisms 
was also analyzed. Due to the lack of genomic information in 
Portunus trituberculatus, the majority of the assembled 
sequences (11.25%) matched genes from microcrustacean 
arthropod Daphnia pulex which have full genomic information. 

Comparison of gene expression among the different 
treatment groups in the current experiment is helpful for 
identification of candidate genes underlying response to salinity 
stress in Portunus trituberculatus. In this study, we detected a 
total of 1,705 unigenes were differentially expressed between 
the comparison of LC vs NC and HC vs NC (qvalue<0.005 & | 
log2 (foldchange)|>1), which significantly more abundant than 
the data of previous study, in which only 417 differentially 
expressed genes were found after two different salt challenges 
(10% and 40%) via cDNA microarray technology [5]. The 
results prove that the next generation sequencing methods 
should be more powerful than cDNA microarray technology in 
expression analysis because it can provide high data output 
and recognize new unigenes or unique isoforms present in 
transcriptome [32]. 

Many fewer differentially expressed unigenes were found in 
the LC Vs. NC (615) compared with that of HC Vs. NC (1,516). 
To further unravel the significantly altered biological processes 
upon salinity stress, the up- and down-regulated genes 
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Figure 5. Eight expression patterns ofthe differentially expressed genes. Cluster I, genes highly up-regulated in low and high 
salinity conditions; Cluster II, genes down-regulated in low and high salinity conditions;Cluster III, genes highly up-regulated in low 
salinity and down-regulated in high salinity challenge; Cluster IV, genes highly up-regulated in high salinity and down-regulated in 
low salinity challenge; Cluster V, genes highly up-regulatedonly in low salinity challenge; Cluster VI, genes down-regulation only in 
low salinity challenge; Cluster VII, genes highly up-regulatedonly in high salinity challenge; Cluster VIII, genes down-regulatedonly 
in high salinity challenge. 

doi: 10.1371/joumal.pone.0082155.g005 



categorized into eight patterns based on expression profiles 
were subjected to the GO term enrichment analysis. The 
results revealed that a total of 454 unigenes were 
enrichmented in 63 processes, among which, the differentially 
expressed unigenes in low salinity stress were enrichmented in 
25 processes (6 processes were induced, 17 processes were 
suppressed and 2 processes were were affected). The 
differentially expressed unigenes in high salinity stress were 
enrichmented in 58 processes (31 processes were induced, 24 
processes were suppressed and 3 processes were were 
affected). Further analysis found, 5 and 38 processes were 
only enrichmented in low salinity stress or high salinity stress 
respectively, another 20 processes were enrichmented in both 
stress. In conclusion, our analysis uncovered: 1, more genes 
and processes were involved in high salinity stress adaptation 
than low salinity stress adaptation. 2, the processed involved in 
low salinity stress were mainly supressed, and the processed 
involved in high salinity stress were mainly induced. 3, hypo- 
osmoregulation and hyper-osmoregulation mechanisms share 
some processes but also has many differences. 

Past research has revealed two major strategies, i.e. the 
"limiting process" and the "compensatory process", that are 
adopted by crustaceans for osmoregulation and both are 
predominately accomplished by the gills [33,34]. A "limiting 
processes a strategy whereby the maintenance of hemolymph 
osmolality/ions is accomplished by adjusting the permeability of 
the boundary structures (e.g. gill mem-branes), which can be a 
highly effective method to reduce ion diffusion and water inux 



rather than solely relying on the more energetically demanding 
mechanisms of ion transport. It has been suggested that the 
mechanisms during long term salinity exposure is often the 
result of gill membrane fatty acid compositional changes [35]. 
In the case, generally membranes containing higher 
proportions of fatty acids with high unsaturation indexes (i.e. 
higher double bonds). In this research, we found at least five 
up-regulation unigenes belong to fatty acid biosynthetic 
process (GO:0006631) in response to salinity stress (Table S4) 
which suggested the strategies of "limiting process" maybe play 
an important role in osmoregulation in Portunus trituberculatus. 

A "compensatory process", on the other hand, is a strategy 
whereby the maintenance of hemolymph osmolality/ions is 
accomplished via the active movement of solutes into or out of 
the hemolymph to counter balance their passive diffusion 
[4,33]. When the environmental osmolality is higher than that of 
the hemolymph, it is termed "hyper-osmotic" and crustaceans 
must compensate for ion influx via hypo-osmoregulation; for 
the opposite, when the environmental osmolality is lower than 
that of the hemolymph, it is termed "hypo-osmotic" and 
crustaceans must compensate for ion loss from their 
hemolymph via hyper-osmoregulation [34]. Portunus 
trituberculatus belongs to weak hyperosmoregulators [36] 
which employs the apical Na + ' K + , 2CI() cotransporter driven by 
the inwardly directed Na + gradient, supplemented by apical K + 
channels, that recycle K + and hyperpolarize the apical 
membrane, creating a negative cell potential that drives Cl~ 
efflux through basal CI channels. Na + uptake may be 
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Figure 6. GO enrichment analysis of the differentially expressed genes. Cluster I, genes highly up-regulated in low and high 
salinity conditions; Cluster II, genes down-regulated in low and high salinity conditions;Cluster III, genes highly up-regulated in low 
salinity and down-regulated in high salinity challenge; Cluster IV, genes highly up-regulated in high salinity and down-regulated in 
low salinity challenge; Cluster V, genes highly up-regulatedonly in low salinity challenge; Cluster VI, genes down-regulation only in 
low salinity challenge; Cluster VII, genes highly up-regulatedonly in high salinity challenge; Cluster VIII, genes down-regulatedonly 
in high salinity challenge. 

doi: 10.1371/joumal.pone.0082155.g006 



augmented through the apical sodium/hydrogen exchanger. 
The resulting outside-positive transepithelial potential drives 
substantial inward paracellular Na + movement across this leaky 
epithelium [36,37]. In our study, ion transport-related genes 
including chloride channel protein, sodium/hydrogen 
exchanger, sodium/glucose cotransporter, and carbonic 
anhydrase were up-regulated in response to low salinity stress. 
In addition, we found that five signifantly up-regulated 
expressed unigenes were enrichmented in anion transport 
process (GO:0006820). The differential expression of this key 
ion transport genes potentially supports this models of decapod 
Crustacea hyper-osmoregulation reported previously [36,37]. 

On the contrary, in hypo-osmoregulation, passive NaCI influx 
is compensated by active NaCI secretion, however, the 
mechanisms of active NaCI secretion in the Crustacea are not 
known[36]. In this studys, we also found ion transport-related 
genes in high salinity stress, but most of them down-regulated 
expression which were different from low salinity stress. These 
down-regulated genes contained Na + , K + , 2CI() cotransporter, 
carbonic anhydrase, calcium-activated chloride channel 
regulator 1, and anion exchange protein etc. Based on 
previous researches [36,37] and our results, we speculated 
that: in hypo-osmoregulation, for maintaining hemolymph 
osmotic balance, Portunus trituberculatus initiatively down- 
reglated the expression of some important ion transport-related 
genes to prevent high concentration of ions of the external 
environment penetrate into the body. However, more detailed 



studies are required to disclose the mechanisms of hypo- 
osmoregulation in Crustacea. 

The aquaporins were considered as important salt 
responsive effectors [38]. The functions of AQP have been 
studied in mammals [39] and fish [40] which suggested that 
epithelial water flux occurs, at least in part, through 
transcellular pathways formed via aquaporins (AQPs), 
however, the functions of aquaporins in osmoregulation was 
little known in Crustacea. Previous researches suggested that 
during acclimation to a hyperosmotic environment, increasing 
renal AQPs expression may play a role in enhancing water 
reabsorption by the kidney [41] and translational gene 
knockdown of AQPs protein reduced water influx in fish [42]. 
Additionally, the water channel activity of AQPs can be 
modulated through phosphorylation and dephosphorylation 
[43-45]. In addition to maintaining ionic homeostasis under salt 
stress conditions, Portunus trituberculatus also need establish 
water homeostasis via aquaporins (AQPs). Our data showed 
that the mRNA expression level of AQPs decreased under 
hypo-osmotic and hyper-osmotic stress conditions. Similar 
results were found in many previous studies under long-term 
saltstress [38,46]. These results may indicated that long-term 
stress induced a reduction of AQPs activities via transcriptional 
level controls and post-translational modifications to protect 
against water currents as well as gill cell swelling and 
shrinkage. However, AQPs activities whether affected by 
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Table 5. GO enrichment analysis of the differentially expressed genes. 



Cluster 


GOID 


Term 


Level of GO 


Number of unigenes 


P-value 


1 


GO:0006030 


chitin metabolic process 


6 


3 


0.001204934 


II 


GO:0006030 


chitin metabolic process 


6 


13 


1 .60E-08 


II 


GO:0051258 


protein polymerization 


7 


6 


0.000898448 


II 


GO:0006813 


potassium ion transport 


8 


4 


0.032189626 


II 


GO:0006979 


response to oxidative stress 


3 


4 


0.000251565 


II 


GO:0006120 


mitochondrial electron transport, NADH to ubiquinone 


9 


3 


0.037277062 


II 


GO:0006875 


cellular metal ion homeostasis 


8 


3 


0.041235929 


II 


GO:0009755 


hormone-mediated signaling pathway 


7 


3 


0.037277062 


II 


GO:0016568 


chromatin modification 


6 


3 


0.043988317 


II 


GO:0071941 


nitrogen cycle metabolic process 


3 


3 


0.005597731 


II 


GO:0006935 


chemotaxis 


4 


2 


0.048209006 


III 


GO:0055114 


oxidation-reduction process 


3 


4 


0.017862676 


III 


GO:0006629 


lipid metabolic process 


3 


3 


0.035596858 


III 


GO:0006546 


glycine catabolic process 


10 


2 


9.13E-06 


III 


GO:0008033 


tRNA processing 


9 


2 


0.002337968 


IV 


GO:0055114 


oxidation-reduction process 


3 


5 


0.043166178 


IV 


GO:0006508 


proteolysis 


5 


3 


0.031323303 


IV 


GO:0002098 


tRNA wobble uridine modification 


12 


2 


0.000102853 


IV 


GO:0000154 


rRNA modification 


10 


2 


0.000190739 


IV 


GO:0009452 


7-methylguanosine RNA capping 


9 


2 


0.000190739 


IV 


GO:0001510 


RNA methylation 


8 


2 


0.000395991 


IV 


GO:0006479 


protein methylation 


8 


2 


0.000770984 


V 


GO:0006820 


anion transport 


6 


5 


5.52E-05 


V 


GO:0006952 


defense response 


3 


5 


0.034822636 


V 


GO:0006614 


SRP-dependent cotranslational protein targeting to membrane 


10 


2 


0.008743629 


VI 


GO:0006508 


proteolysis 


5 


8 


0.003397869 


VI 


GO:0006030 


chitin metabolic process 


6 


3 


0.01432079 


VI 


GO:0006888 


ER to Golgi vesicle-mediated transport 


7 


2 


0.037838935 


VI 


GO:0009878 


nodule morphogenesis 


6 


2 


0.021455881 


VI 


GO:0001539 


ciliary or flagellar motility 


5 


2 


0.001800326 


VII 


GO:0055114 


oxidation-reduction process 


3 




0.000952872 


VII 


GO:0006464 


cellular protein modification process 


6 


50 


0.000230686 


VII 


GO:0006508 


proteolysis 


5 


40 


0.000397024 


VII 


GO:0016310 


phosphorylation 


5 


31 


0.001296757 


VII 


GO:0008610 


lipid biosynthetic process 


4 


19 


0.003426418 


VII 


GO:0009966 


regulation of signal transduction 


5 


17 


0.041422297 


VII 


GO:0007264 


small GTPase mediated signal transduction 


6 




0.000113889 


VII 


GO:0006030 


chitin metabolic process 


6 


11 


0.006871104 


VII 


GO:0009108 


coenzyme biosynthetic process 


5 


10 


0.001612753 


VII 


GO:0006888 


ER to Golgi vesicle-mediated transport 


7 


7 


0.021852563 


VII 


GO:0006006 


glucose metabolic process 


7 


6 


0.038801546 


VII 


GO:0042398 


cellular modified amino acid biosynthetic process 


9 


4 


0.04770803 


VII 


GO:0035023 


regulation of Rho protein signal transduction 


8 


4 


0.008616499 


VII 


GO:0032446 


protein modification by small protein conjugation 


8 


4 


0.039101467 


VII 


GO:0007339 


binding of sperm to zona pellucida 


7 


4 


0.000990193 


VII 


GO:0006099 


tricarboxylic acid cycle 


7 


3 


0.047820393 


VII 


GO:0006000 


fructose metabolic process 


7 


3 


0.007238233 


VII 


GO:0001503 


ossification 


3 


3 


0.043900004 


VII 


GO:0009089 


lysine biosynthetic process via diaminopimelate 


11 


2 


0.04770803 


VII 


GO:0015813 


L-glutamate transport 


11 


2 


0.018737173 


VII 


GO:0006425 


glutaminyl-tRNA aminoacylation 


11 


2 


0.018737173 


VII 


GO:0009168 


purine ribonucleoside monophosphate biosynthetic process 


10 


2 


0.041097287 


VII 


GO:0034968 


histone lysine methylation 


10 


2 


0.03486768 


VII 


GO:0046058 


cAMP metabolic process 


10 


2 


0.02365927 
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Table 5 (continued). 



Cluster 


GOID 


Term 


Level of GO 


Number of unigenes 


P-value 


VII 


GO:0006570 


tyrosine metabolic process 


9 


2 


0.014309718 


VII 


GO:0006541 


glutamine metabolic process 


9 


2 


0.04770803 


VII 


GO:0006784 


heme a biosynthetic process 


8 


2 


0.010408479 


VII 


GO:0009231 


riboflavin biosynthetic process 


7 


2 


0.029045762 


VII 


GO:0030534 


adult behavior 


4 


2 


0.03486768 


VIII 


GO:0055114 


oxidation-reduction process 


3 


27 


0.001223294 


VIII 


GO:0006508 


proteolysis 


5 


20 


1.19E-05 


VIII 


GO:0005975 


carbohydrate metabolic process 


3 


17 


6.23E-05 


VIII 


GO:0006030 


chitin metabolic process 


6 


10 


2.79E-06 


VIII 


GO:0006544 


glycine metabolic process 


9 


4 


0.000286422 


VIII 


GO:0006631 


fatty acid metabolic process 


8 


4 


0.013582145 


VIII 


GO:0042742 


defense response to bacterium 


5 


4 


0.020716006 


VIII 


GO:0006825 


copper ion transport 


9 


3 


0.002850419 


VIII 


GO:0007218 


neuropeptide signaling pathway 


7 


3 


0.042788588 


VIII 


GO:0006979 


response to oxidative stress 


3 


3 


0.002600038 


VIII 


GO:0006334 


nucleosome assembly 


10 


2 


0.01916756 


VIII 


GO:0009396 


folic acid-containing compound biosynthetic process 


10 


2 


0.01916756 


VIII 


GO:0006563 


L-serine metabolic process 


9 


2 


0.00077586 


VIII 


GO:0009084 


glutamine family amino acid biosynthetic process 


9 


2 


0.005774778 


VIII 


GO:0008272 


sulfate transport 


8 


2 


0.011002118 


VIII 


GO:0007631 


feeding behavior 


3 


2 


0.020714011 



doi: 10.1371/joumal.pone.0082155.t005 



phosphorylation in Portunus trituberculatus should be further 
verified. 

It has been reported that Na + ,K + -ATPase is the driving force 
in establishing an ion gradient across the epithelial cell 
membrane in marine crabs [47], In this study, Na + ,K + -ATPase 
P-subunit showed significant down-regulation in high and low 
salinity conditions, which suggested that Na + ,K + -ATPase might 
play a very important role in salinity adaptation. Previous 
studies have shown that Na + ,K + -ATPase p-subunit showed 
significant up-regulation during low salinity challenge and 
significant down-regulation during high salinity challenges [5], 
which was different from our data. This may be due to different 
salinity challenged in the two studies. Although lots of literature 
showed that gene expression levels of Na + ,K + -ATPase a 
subunit were highly up-regulated during salinity stress [48-50], 
significant up-/down-regulation of Na + ,K + -ATPase a subunit 
was not detected in our transcriptome data, and similar results 
were also found in Xu's studies [5]. None of the significant 
gene expression changes of Na + ,K + -ATPase a subunit during 
salinity challenges do not negate its role for osmoregulation in 
swimming crab, and more future researches should be done 
toward Na + ,K + -ATPase a subunit. 

Previous studies have shown that the gill apical V-H + - 
ATPase in freshwater crabs was shown to be involved in ion 
regulation, and marine ones generally showed a cytoplasmic V- 
H + -ATPase distribution, which did not participate in 
osmoregulation [51]. Similar suggestions have already been 
proposed in studies of the gills of Eriocheir sinensis and ilea 
tangeri [52-54]. Portunus trituberculatus belong to typical 
euryhaline marine crab species, however, it is interesting that 
four V-H+-ATPase subunits were showed significantly different 



expresstion in high salinity stress, which suggested V-H + - 
ATPase might play a very important role during hypo- 
osmoregulation. Subcellular localization and function of V-H + - 
ATPase deserve further study in Portunus trituberculatus. 

The present investigation suggests that the degradation of 
hemolymph or muscle protein to FAAs or de novo synthesis of 
free amino acids (FAA) may serve as a mechanism to 
compensate or alleviate the effects of ion influx under salinity 
stress [55], and that glycine, proline and alanine widely function 
as osmoeffectors in a number of crustacean species [56-59]. It 
must be noted, however, such a response is not universal to all 
crustaceans [60]. In this study, many differentially expressed 
unigenes were enrichmented in amino acid metabolism and 
synthesis processes. Among which, glycine catabolic process 
(GO:0006546) was induced in low salinity stress, and in high 
salinity stress, glycine metabolic process (GO:0006544) and L- 
serine metabolic process (GO:0006563) were suppressed, but 
lysine biosynthetic process (GO:0009089) was induced. 
Interestingly, we aslo found many differentially expressed 
unigenes enriched in proteolysis process (GO:0006508), which 
down-regulated (1 1 genes) in low salinity stress and mainly up- 
regulated (43 up-regulated genes and 20 down-regulated 
genes) in high salinity stress. This results uncovered that, in 
high salinity stress, to balance the osmotic pressure inside and 
outside the body, Portunus trituberculatus increase levels of 
FAAs in vivo via the way of increasing protein degradation, 
decreaseing amino acid catabolism, and increasing synthesis 
of certain non-essential amino acids, and in low salinity stress, 
Portunus trituberculatus takes the opposite approach to 
hypertonic regulation. 
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In addition to the known processes associated with 
osmoregulation in Crustacea., there were some other 
processes enriched which were little known in osmoregulation. 
The top three processes based on the number of unigenes 
were oxidation-reduction process (GO:00551 14), cellular 
protein modification process (GO:0006464) and 
phosphorylation (GO:0016310), and the top three processes 
based on the enrichment level (p-value, excluding the 
processes of less than three unigenes) were chitin metabolic 
process (GO:0006030), carbohydrate metabolic process (GO: 
0005975) and small GTPase mediated signal transduction 
(GO:0007264). 

The oxidation-reduction process is essential in maintaining 
cellular homeostasis. Under physiologic conditions, cells 
maintain oxidation-reduction balance through generation and 
elimination of reactive oxygen/nitrogen species (ROS/RNS) 
[61]. Normally, the oxidation-reduction process homeostasis 
ensures that the cells respond properly to endogenous and 
exogenous stimuli. However, when the oxidation-reduction 
process homeostasis is disturbed, oxidative stress may lead to 
aberrant cell death and contribute to disease development [61]. 
Both exogenous and endogenous sources contribute to the 
formation of intracellular ROS/RNS. This reseach suggested 
that exposure to salinity stress, especially the high salinity 
stress has been shown to active the oxidation-reduction 
process in Portunus trituberculatus in order to balance the level 
of ROS/RNS. Of course, there was another possibility that the 
oxidation-reduction process in Portunus trituberculatus has 
been disturbed in this high salinity stress. Both cellular protein 
modification and phosphorylation process belong to the scope 
of post-translational modifications. 

Posttranslational modifications modulate the activity of most 
eukaryote proteins [62] which, manifest as chemical 
modifications that occur on amino acid side chains in a site- 
specific manner. They can temporarily or permanently change 
the fate of the protein by enhancing the functionality and/or 
stability of the target protein through the recruitment of auxiliary 
factors, change the proteins' cellular localisation or signal the 
most terminal fate, proteasomal degradation [63]. Our data 
showed that many gene were enrichmented in cellular protein 
modification (50 unigenes ) and phosphorylation process (31 
unigenes ), suggested that posttranslational modifications play 
an important role in osmoregulation of Portunus trituberculatus, 
and further research is worthwhile. 

Many differentially expressed unigenes (38) enriched in 
chitin metabolic process with the maximum enrichment level (p- 
value: 1.60E-08) were aslo found in our data, which suggested 
chitin metabolic process may be involved in osmoregulation or 
affected by salinity stress. Further analysis found there were 4 
chitinase genes. In crustaceans, during the molting cycle, 
chitinase dissolves chitin in the old exoskeleton into more 
soluble forms, which can then be partially reabsorbed into the 
body and used to synthesize the new exoskeleton [64]. 
Therefore, the chitinase are indispensable for crustaceans. By 
now, many genes encoding chitinases have been isolated and 
characterized in crustaceans [64-67] due to their importance in 
growth, development and immunity. However, chitinase with 
function of osmoregulation has not been reported in Crustacea. 



More detailed research can be done according to the 
connections between the salinity changes and chitin metabolic 
process in the future. 

Small GTP-binding genes play a crucial regulatory role in a 
number of cellular processes in both plants and ani-mals, such 
as vesicle-mediated intracellular trafficking, signal transduction, 
cytoskeletal organization, and cell division [68]. Differential 
expression of small GTP-binding proteins has been previously 
reported in various abiotic stresses in plant [68], animal [69], 
and bacteria [70], however, the functions of Small GTP-binding 
genes was little known in Crustacea. 15 unigenes were 
enriched in GTPase mediated signal transduction process with 
a relatively low p-value implied that GTPase mediated signal 
transduction process plays an important role in Portunus 
trituberculatus salinity adaptation. In addition, 17 up-regulated 
unigenes were enriched in carbohydrate metabolic process in 
high salinity stress suggested that hypo-osmoregulation of 
Portunus trituberculatus depends on energy consumption. This 
result have also been reported in Xu's studies [5] and further 
experiments are needed to verify the hypothesis. 

It should be noted that, the processes of response to 
oxidative stress (GO:0006979), defense response (GO: 
0006952), response to oxidative stress (GO:0006979) and 
defense response to bacterium (GO:0042742) were enriched in 
salinity stess, and lots of stress-related and immunity-related 
genes such as HSP, cathepsin, lectin, serine protease, and 
peroxiredoxin show significantly up-regulated or down- 
regulation expression responding to low or high salinity stress, 
which suggested that salinity dilution or elevation appeared to 
invoke a classic "stress response" and "immunity response" at 
the transcriptional level in Portunus trituberculatus. This results 
similar to Xu et al.(2011)'s results [5] and different from Towle 
et al.(2011)'s results [71], which might reflect that Portunus 
trituberculatus possesses lower salinity adaptability than C. 
maenas, and it also suggested that the two crab species might 
have different signal pathways against the salinity "stress"[5]. 

It is noteworthy that the present study identified a large 
number of transcripts (526 unigenes) significantly upregulated 
or downregulate during salinity stress for which no annotation 
was readily available. This new genes provide a wealth of 
reference data to further research the mechanism of 
osmoregulation in Crustacea. 

Conclusions 

In general, our work represents the first report of the 
utilization of the next generation sequencing techniques for the 
study of osmoregulation in Portunus trituberculatus. Through 
more than 18.16G clean bases obtained, we established 
transcriptome data set for Portunus trituberculatus subjected to 
salinity stress, and the data we generated could enrich on 
genomic resources of this non-model organism. Based on the 
comparison of gene expression, a substantial number of genes 
were found to be modified by salinity stress which 
demonstrated the complexity of salinity adaptation mechanism 
in the crab. Our results revealed a few important salinity 
acclimation pathways via enrichment analysis, which may be 
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helpful in understanding the molecular basis of osmoregulation 
in Portunus trituberculatus. 
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Table S1. Primer sequences and product size of target 
and reference genes used for quantitative RT-PCR. 

(XLSX) 

Table S2. List of differently expressed genes from 
Portunus trituberculatus in response to low salinity 
challenged (LC) versus non-challenged (NC). 
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Table S3. List of differently expressed genes from 
Portunus trituberculatus in response to high salinity 
challenged (HC) versus non-challenged (NC). 
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Figure S1. COG Classification of the unigenes. Possible 
functions of 11528 unigenes were classified and subdivided 
into 26 COG categories. 
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Figure S2. GO classification of all unigenes. Most unigenes 
canbe divided into three major categories, including biological 
process, cellular component, and molecular function. 
(JPG) 

Figure S3. KEGG Classification of the unigenes. 5419 
unigenes were assigned into 31 KEGG pathways.A, Cellular 
Processes; B, Enviromental Information Processing; C, Genetic 
Information Processing; D, Metabolism; E, Organismal 
Systems. 
(PNG) 
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